function [error,error_25]=model5d_driftcorrection_toolbox(patient,bellows_volt_drifted, flow_drifted, ii,jj,kk);
    
    %| 5D model code, for use in the drift correction segment of the
    %workflow
    %------------------------------------------------------------------------
    %|      Dependancies;                                                     |
    %|           read_dvf_elastix_slice_toolbox                   |
    %------------------------------------------------------------------------
    %   This file is part of the
    %   5D-Novel4DCT Toolbox  ("Novel4DCT-Toolbox")
    %   DH Thomas, Ph.D
    %   University of California, Los Angeles
    %   Contact: mailto:dhthomas@mednet.ucla.edu
    %------------------------------------------------------------------------
    % $Author: DHThomas $	$Date: 2014/04/01 10:23:59 $	$Revision: 0.1 $

    indexs_vec=repmat(kk,size(jj,2),1);
    indexs_vec=indexs_vec(:)';
    
    X1=bellows_volt_drifted(indexs_vec,patient.run_scans)';
    X2=flow_drifted(indexs_vec,patient.run_scans)';
    
    [xgrid, ygrid]=meshgrid(1:size(X1,2),1:size(patient.run_scans,2));
    
    error=zeros(patient.dim(1),size(jj,2),size(kk,2));
    error_25=zeros(patient.dim(1),size(jj,2),size(kk,2));
    
    
    for slice=ii;%round(dim(1)/2);%1:dim(1);%par
        
        [ dvf]=read_dvf_elastix_slice_toolbox(patient.dvf_folder,slice);
        
        % subsample the DVF slices, by [run_scans, jj,kk];
        dvf_u_vec=single(dvf.u(patient.run_scans,jj,kk));dvf_u_vec=dvf_u_vec(:,:);
        dvf_v_vec=single(dvf.v(patient.run_scans,jj,kk));dvf_v_vec=dvf_v_vec(:,:);
        dvf_w_vec=single(dvf.w(patient.run_scans,jj,kk));dvf_w_vec=dvf_w_vec(:,:);
        
        tempfit_variables_x=zeros(3,length(X1));%length(initparams));
        tempfit_variables_y=zeros(3,length(X1));%length(initparams));
        tempfit_variables_z=zeros(3,length(X1));%length(initparams));
        
        yfit_x=zeros(size(patient.run_scans,2),length(X1));%length(initparams));
        yfit_y=zeros(size(patient.run_scans,2),length(X1));%length(initparams));
        yfit_z=zeros(size(patient.run_scans,2),length(X1));%length(initparams));
        
        constant=zeros(length(jj),length(kk),3);
        alpha=zeros(length(jj),length(kk),3);
        beta=zeros(length(jj),length(kk),3);
        %     gamma=zeros(length(jj),length(kk),3);
        
        for vox=1:length(X1);%%loop #1 all 25 scans
            
            E_vox=[ones(size(X1,1),1) X1(:,vox) X2(:,vox)];% X3(:,vox)];
            
            tempfit_variables_x(:,vox)=E_vox\dvf_u_vec(:,vox);
            tempfit_variables_y(:,vox)=E_vox\dvf_v_vec(:,vox);
            tempfit_variables_z(:,vox)=E_vox\dvf_w_vec(:,vox);
            
            yfit_x(:,vox)=tempfit_variables_x(:,vox)'*E_vox';
            yfit_y(:,vox)=tempfit_variables_y(:,vox)'*E_vox';
            yfit_z(:,vox)=tempfit_variables_z(:,vox)'*E_vox';
            
        end
        
        error_u = bsxfun(@minus,yfit_x,dvf_u_vec);
        error_v = bsxfun(@minus,yfit_y,dvf_v_vec);
        error_w = bsxfun(@minus,yfit_z,dvf_w_vec);
        
        %     error_tmp=mean(bsxfun(@hypot,error_w,bsxfun(@hypot,error_u,error_v)));
        error_vec=bsxfun(@hypot,error_w,bsxfun(@hypot,error_u,error_v));
        error_tmp_25=mean(error_vec);
%         [c scan_index]=sort(error_vec);
%         %     leave_out=scan_index(21:25,:);
%         leave_in=scan_index(1:size(patient.run_scans,2),:);
%         leave_in_index=sub2ind(size(X1),leave_in,xgrid);
%         X1_2=X1(leave_in_index);% X1_2=reshape(X1_2,20,size(X1,2));%clear X1;
%         X2_2=X2(leave_in_index);% X2_2=reshape(X2_2,20,size(X1,2));%clear X2;
%         dvf_u_vec_2=dvf_u_vec(leave_in_index); %dvf_u_vec_2=reshape(dvf_u_vec_2,20,size(X1,2));%clear X1;
%         dvf_v_vec_2=dvf_v_vec(leave_in_index); %dvf_v_vec_2=reshape(dvf_v_vec_2,20,size(X1,2));%clear X1;
%         dvf_w_vec_2=dvf_w_vec(leave_in_index); %dvf_w_vec_2=reshape(dvf_w_vec_2,20,size(X1,2));%clear X1;
%         
%         tempfit_variables_x=zeros(3,size(X1_2,2));%length(initparams));
%         tempfit_variables_y=zeros(3,size(X1_2,2));%length(initparams));
%         tempfit_variables_z=zeros(3,size(X1_2,2));%length(initparams));
%         
%         yfit_x=zeros(size(leave_in,1),size(X1_2,2));%length(initparams));
%         yfit_y=zeros(size(leave_in,1),size(X1_2,2));%length(initparams));
%         yfit_z=zeros(size(leave_in,1),size(X1_2,2));%length(initparams));
%         
%         
%         for vox=1:size(X1_2,2);%%loop #2 minus the 5 scans with largest errors
%             
%             E_vox=[ones(size(X1_2,1),1) X1_2(:,vox) X2_2(:,vox)];% X3(:,vox)];
%             
%             tempfit_variables_x(:,vox)=E_vox\dvf_u_vec_2(:,vox);
%             tempfit_variables_y(:,vox)=E_vox\dvf_v_vec_2(:,vox);
%             tempfit_variables_z(:,vox)=E_vox\dvf_w_vec_2(:,vox);
%             
%             yfit_x(:,vox)=tempfit_variables_x(:,vox)'*E_vox';
%             yfit_y(:,vox)=tempfit_variables_y(:,vox)'*E_vox';
%             yfit_z(:,vox)=tempfit_variables_z(:,vox)'*E_vox';
%             
%         end
%         
%         error_u = bsxfun(@minus,yfit_x,dvf_u_vec_2);
%         error_v = bsxfun(@minus,yfit_y,dvf_v_vec_2);
%         error_w = bsxfun(@minus,yfit_z,dvf_w_vec_2);
%         
%         %     error_tmp=mean(bsxfun(@hypot,error_w,bsxfun(@hypot,error_u,error_v)));
%         error_vec=bsxfun(@hypot,error_w,bsxfun(@hypot,error_u,error_v));
%         error_tmp=nanmean(error_vec);
%         error_tmp=sum(error_vec(1:patient.scans,:),1);
%         
%         error_25(slice,:,:)=reshape(error_tmp_25,size(jj,2),size(kk,2));
        error(slice,:,:)=reshape(error_tmp_25,size(jj,2),size(kk,2));
        
        %     error_index(jj,kk,:)=reshape(scan_index())
        
%         constant(jj,kk,1)=reshape(tempfit_variables_x(1,:),size(jj,2),size(kk,2));
%         constant(jj,kk,2)=reshape(tempfit_variables_y(1,:),size(jj,2),size(kk,2));
%         constant(jj,kk,3)=reshape(tempfit_variables_z(1,:),size(jj,2),size(kk,2));
%         %
%         alpha(jj,kk,1)=reshape(tempfit_variables_x(2,:),size(jj,2),size(kk,2));
%         alpha(jj,kk,2)=reshape(tempfit_variables_y(2,:),size(jj,2),size(kk,2));
%         alpha(jj,kk,3)=reshape(tempfit_variables_z(2,:),size(jj,2),size(kk,2));
%         
%         beta(jj,kk,1)=reshape(tempfit_variables_x(3,:),size(jj,2),size(kk,2));
%         beta(jj,kk,2)=reshape(tempfit_variables_y(3,:),size(jj,2),size(kk,2));
%         beta(jj,kk,3)=reshape(tempfit_variables_z(3,:),size(jj,2),size(kk,2));
        
        %     gamma(jj,kk,1)=reshape(tempfit_variables_x(4,:),size(jj,2),size(kk,2));
        %     gamma(jj,kk,2)=reshape(tempfit_variables_y(4,:),size(jj,2),size(kk,2));
        %     gamma(jj,kk,3)=reshape(tempfit_variables_z(4,:),size(jj,2),size(kk,2));
        
%         if save_files==1
%             cd(folder)
%             cd(folder_output)
%             iSaveX( sprintf( 'constant_%d', slice ), constant );
%             iSaveX( sprintf( 'alpha_%d', slice ), alpha );
%             iSaveX( sprintf( 'beta_%d', slice ), beta );
%             %         iSaveX( sprintf( 'gamma_%d', slice ), gamma );
%         end

        % toc
        %     if size(ii,2)>1
        %     ppm.increment();
        %     end
    end
    